/*
 *  Matrix.h
 *  NSE_SOLVER
 *
 *  Created by Patrick Zulian on 12/28/10.
 *  Copyright 2010 ICS - USI - Lugano. All rights reserved.
 *
 */

#ifndef PIDAMR_MATRIX_H
#define PIDAMR_MATRIX_H

#include <iostream>

#include "Numerics.h"
#include "Array.h"
#include "Vector.h"
#include "MatrixHelper.h"
#include "Range.h"
#include "MatrixHandle.h"

#include "AccessInterfaceHandles.h"

namespace nsesolver {
	
	template< typename T = double, class Storage = Array< T >, class MImpl = MatrixHelper< T >  >
	class Matrix {
	private:		
		int nX, nY;
		Storage data;
		
	public:
		typedef T Scalar;
		
		Matrix(int nX, int nY)
		: nX(nX), nY(nY), data(nX*nY)
		{
			
		}
        
		
		Matrix(const Matrix &other) 
		: nX(other.nX), nY(other.nY), data(other.data)
		{
			
		}
		
		Matrix & operator=(const Matrix &other)
		{
			if(this == &other)
				return *this;
			
			nX = other.nX;
			nY = other.nY;
			data = other.data;
			return *this;
		}
        
		//Output
		friend
		std::ostream & operator <<(std::ostream &os, const Matrix &matrix)
		{
			
            //os << i << ".\t";
            for (int j = 0; j < matrix.nY; j++) {
                for (int i = 0; i < matrix.nX; i++) {
					if (Math::abs< T >(matrix(i, j)) < 10e-64) {
						os << 0 << " ";
						continue;
					}
					
					os << matrix(i,j) << " ";
					
				}
                
				os << std::endl;
                
			}
			return os;
		}
		
		
		//Access Functions
		inline
		T & operator()(int x, int y)
		{
			NSE_SOLVER_ASSERT( x < nX );
			NSE_SOLVER_ASSERT( x >= 0 );
			NSE_SOLVER_ASSERT( y < nY );
			NSE_SOLVER_ASSERT( y >= 0 );
			
			
			return data[ x+ y * xmax() ];
		}
		
		inline
		const T & operator() (int x, int y) const 
		{
            NSE_SOLVER_ASSERT( x < nX );
			NSE_SOLVER_ASSERT( x >= 0 );
			NSE_SOLVER_ASSERT( y < nY );
			NSE_SOLVER_ASSERT( y >= 0 );
			
			return data[ x  + y * xmax() ];
		}
		
		inline
		const T & operator[] (int index) const 
		{
			return data[index];
		}
		
		inline
		T & operator[] (int index)  
		{
			return data[index];
		}
		
		inline
		const T & operator[] (const Limit &limit) const 
		{
			NSE_SOLVER_UNUSED(limit);
			return data[size() -1];
		}
		
		inline
		T & operator[] (const Limit &limit)  
		{
			NSE_SOLVER_UNUSED(limit);
			return data[size() -1];
		}
		
		
		
		
		inline
		int ymax() const
		{
			return nY;
		}
		
		inline
		int xmax() const
		{
			return nX;
		}
		
		inline
		int size() const
		{
			return xmax() * ymax();
		}
        
		
		const MatrixHandle< Matrix > operator()(const Range &xRange, const Range &yRange) 
		{
			return MatrixHandle< Matrix >(*this, xRange, yRange);
            
		}
        
		Matrix & operator = (const MatrixHandle< Matrix > &handle)
		{
			return *this = handle; 
			
		}
		
		Matrix & allSet(const T value)
		{
			data.allSet(value);
			return *this;
		}
		
		
		
		//Test functions
		bool equals(const Matrix &other) const
		{
			return MImpl::equals(*this, other);
		}
		
		bool operator==(const Matrix &other) const
		{
			return MImpl::equals(*this, other);
		}
        
		//Matrix-scalar operations
		Matrix & operator *= (T scaleFactor)
		{
			MImpl::inPlaceScale(*this, scaleFactor);
			return *this;
		}
        
		
		friend
		const Matrix operator * (const Matrix &matrix, T scaleFactor)  
		{
			return Matrix(matrix) *= scaleFactor;
		}
		
		friend const Matrix operator * (T scaleFactor, const Matrix & matrix) 
		{ 
			return matrix * scaleFactor;
		}
		
        
		inline
		const Matrix operator -() const 
		{
			return (*this) * static_cast< T > (-1);
			
		}
		
		
		
		
		const T getMax()
		{
			using std::max;
			
			T result = (*this)[0];
			for (int i = 1; i < size(); i++) {
				result = max(result, (*this)[i]);
			}
			return result;
		}
		
		
		
		const T getMin()
		{
			using std::min;
			
			T result = (*this)[0];
			for (int i = 1; i < size(); i++) {
				result = min(result, (*this)[i]);
			}
			return result;
		}
		
		
		const T getAbsMax() const
		{
			using std::max;
			
			T result = Math::abs( (*this)[0] );
			for (int i = 1; i < size(); i++) {
				result = max(result, Math::abs( (*this)[i]) );
			}
			return result;
		}
		
		
		
		const T getAbsMin() const
		{
			using std::min;
			
			T result = Math::abs( (*this)[0] );
			for (int i = 1; i < size(); i++) {
				result = min(result, Math::abs( (*this)[i]) );
			}
			return result;
		}
		
        
	};
	
	
	
}

#endif
